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Resumen 


Se aplica el método de los volúmenes finitos y se generaliza la estrategia ultimate para 
la solución de la ecuación advección-difusión-reacción con parámetros variables, a 
fin de obtener resultados con difusión numérica despreciable. El modelo se verifica 
mediante el cálculo de problemas con solución analítica conocida, mostrando su 
aptitud para simular el transporte difusivo y advectivo de propiedades escalares en 
el agua. Se aplica para calcular la temperatura en un tramo de 32 km del río Itata 
bajo, ubicado en Chile Central, donde se midió la temperatura del río y las variables 
que controlan el régimen térmico en varios puntos. Se presentan resultados para 
cuatro periodos de doce días de duración, correspondientes a escenarios de invierno, 
primavera, verano y otoño. La temperatura media medida en los escenarios de 
invierno y verano es constante, igual a 9.6 y 24.0 *C, respectivamente, y no cambia 
a lo largo del tramo. En los escenarios medidos de otoño y primavera se observa 
un leve aumento a lo largo del tramo igual a 0.45 y 0.74 *C, respectivamente. La 
amplitud de la oscilación térmica diaria medida en el escenario de invierno es 
despreciable, igual a 0.7 *C y muy significativa en el de verano, igual a 8.4 *C. Se 
observa que sólo en el escenario medido de verano dicha amplitud aumenta a lo 
largo del tramo de estudio, a razón de 0.175 “C /km. El modelo desarrollado predice 
de manera adecuada el complejo comportamiento que exhibe la temperatura 
medida durante los cuatro escenarios analizados, constituyendo una herramienta 
de simulación con capacidad predictiva. 


Palabras clave: régimen térmico, métodos numéricos, esquemas de discretización, 
río Itata, Chile. 


Introducción 


El régimen térmico de un río controla pro- 
piedades agua, 
densidad y viscosidad, afectando la turbulencia, 


fundamentales del como 
y con ello la capacidad de transporte y mezcla. 
Es relevante para el ecosistema acuático y la 
calidad del agua, ya que condiciona la presencia 
y distribución de especies, la concentración de 
saturación de gases disueltos y la velocidad de 
reacciones químicas. Caissie (2006) presentó 
una revisión del régimen térmico y sus alcances 
para el ecosistema fluvial. Las represas y 
descargas de efluentes térmicos, por ejemplo, 
provenientes de aguas de enfriamiento de cen- 


trales termoeléctricas y nucleares, pueden alterar 
significativamente el régimen térmico (ver e.g., 
Val et al., 2006). 

Al escurrir, el agua de un río cambia su 
temperatura debido a un intercambio de calor 
con la atmósfera en su superficie, y con el 
lecho en su cauce. Las fuentes y/o sumideros 
de calor son la radiación solar, radiación 
atmosférica, radiación emitida por el agua, 
calor por evaporación, calor convectivo, y calor 
de la lluvia y del lecho (EPA, 1985; Thomann y 
Mueller, 1987; Kim y Chapra, 1997). Alexander 
y Caissie (2003) mostraron que el intercambio 
de calor con la atmósfera explica en más de un 
85% la temperatura del agua, mientras que el 
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intercambio de calor con el lecho explica del 
orden de un 10%, siendo despreciables los 
aportes de calor por lluvia, el intercambio de 
calor con aguas subterráneas y el efecto de la 
sombra producida por la vegetación ribereña, 
especialmente en ríos relativamente grandes. 
Tung et al. (2007) propusieron un modelo 
matemático basado en la ley de Beer para 
incorporar el efecto de la sombra provocada 
por la vegetación ribereña sobre la temperatura 
del río y lo aplicaron al estero GaoShan ubicado 
en Taiwán. La temperatura del agua suele ser 
sensible a la vegetación ribereña en ríos de bajo 
orden, donde el ancho del área sombreada es 
significativo respecto al ancho del río. 

Las fuentes y los sumideros de calor at- 
mosférico pueden determinarse a partir de 
mediciones meteorológicas de cuatro variables 
fundamentales: 1) radiación solar incidente, 2) 
temperatura ambiente, 3) velocidad del viento, 
y 4) humedad relativa del aire. Sinokrot y Stefan 
(1994), y Evans et al. (1998) mostraron que la 
radiación solar incidente y la temperatura 
ambiente son las variables que más afectan a la 
temperatura del agua. 

La distribución espacio-temporal de la 
temperatura de un río no controlado presenta 
cambios de difícil predicción, gobernados 
por ecuaciones no lineales, cuya solución 
en muchos casos requiere de esquemas de 
discretización libres de difusión numérica 
significativa. Kim y Chapra (1997) aplicaron 
un esquema implícito de diferencias finitas 
a la solución de la ecuación del transporte 
de calor en forma no conservativa, mientras 
que Younus et al. (2000) aplicaron el esquema 
predictor-corrector de MacCormack. Westhoff 
et al. (2007) evaluaron la calidad de un modelo 
numérico para el cálculo de temperatura de un 
tramo de 600 m, comparando sus resultados 
con mediciones de alta resolución espacial y 
temporal obtenidas mediante el uso de fibra 
Óptica. Link (2006) presentó un experimento 
donde 
numérica generada por varios esquemas de 


numérico, comparó la difusión 


discretización como diferencias centradas, 


upwind, híbrido y quick. Adicionalmente a 
estos esquemas, pueden seguirse estrategias 
que, en función de la solución local, deciden 
la aplicación de esquemas de diferente orden 
en forma ponderada para limitar los flujos. Por 
ejemplo, Gaskell y Lau (1988) propusieron la 
estrategia smart; Leonard (1991), la estrategia 
ultimate; Woodfield et al. (2004), la estrategia 
wsn; y Qian et al. (2007), la estrategia phib. 

En el presente artículo se presenta un 
modelo conservativo en volúmenes finitos para 
el cálculo de la distribución espacio-temporal 
de la temperatura de un río. Se extiende la 
estrategia ultimate al caso en que la ecuación 
gobernante tiene coeficientes variables. El 
modelo desarrollado se verifica mediante 
problemas con solución analítica conocida y se 
aplica para el cálculo del régimen térmico del 
río Itata bajo durante cuatro periodos de doce 
días de duración cada uno, correspondientes 
a escenarios de invierno, primavera, verano y 
otoño. 


Ecuaciones gobernantes 
La distribución espacio-temporal de la 
temperatura en un río bien mezclado en 
sus direcciones transversal y vertical queda 


descrita por la ecuación advección-difusión- 
reacción: 


AT, (QT) a 
ot OX Ox 


donde A [m?] es el área de la sección trans- 
versal; T [*C], la temperatura media en la sec- 
ción transversal; t [s], el tiempo; Q [m's”], 
el caudal; x [m], la distancia en sentido lon- 
gitudinal; D, [m*s”], el coeficiente de dispersión 
longitudinal; H, [W/m?], el flujo de energía 
térmica a través de la superficie libre y del 
lecho; W [m], el ancho del canal; C, [J kg* *C”], 
el calor específico del agua, y p,, [kg m”] es la 
densidad del agua. El coeficiente de dispersión 
longitudinal, D, [m?s] se calculó mediante la 
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fórmula válida para ríos de pendiente escasa 
propuesta por Vargas y Ayala (2001), pudiendo 
variar en el espacio y en el tiempo: 


w -1.86 2 
e, pa A 
Ry U,RE 


R,, [m] es el radio hidráulico y U. [m s”] es la 
velocidad de corte. El flujo superficial neto de 
energía térmica H,, despreciando aportes por 
lluvia, lecho y sombra, corresponde a: 


Q) 


HA, = Hyy + H y H - HH (3) 
donde H,,, es la radiación neta de onda corta; 
Hyw la radiación atmosférica neta de onda 
larga; H,, la radiación de onda larga emitida 
por el agua; H,, la tasa de transferencia de calor 
por conducción; yH,, la tasa de calor producto 
de la evaporación. 


Radiación neta de onda corta, H,,, 


La radiación neta de onda corta corresponde a: 


SW su Ñ tc (4) 


incidente 


La radiación solar incidente, H , se 
Ñ A Ñ ES indicente 

mide directamente mediante un radiómetro. La 
radiación solar reflejada, H , se calculó 


SWreflejada 
como: 


sw; 


(5) 


donde R, es la fracción de la radiación solar 
que es reflejada en la superficie del agua: 


b 
180 ] (6) 


Ema La 
T 


donde a es la altitud solar en radianes, y a y b 
son coeficientes que dependen de la cobertura 
nubosa. De acuerdo con Martin y McCutcheon 
(1999): 


Si Cb=>=0.9 
Si 0.5<Cb<0.9 
Si 0.1<Cb<0.5 
Si Cb<0.1 


>4=0.33 b=-0.45 
=>4=0.95 b=-0.75 
=>4=220 b=-0.97 
=4=1.18 b=-0.77 


(7) 


La cobertura nubosa se determinó como: 


Ch=1- Hrs sol... 
Hrs sol 


teóricas 


(8) 


Las horas de sol reales se determinaron a 
partir de la radiación solar incidente medida, 
considerando que hay sol cuando durante la 
hora se exceden los 120 W/m?. Cb se determinó 
para cada día. 


Radiación atmosférica neta de onda larga, 
Hwy 

La radiación atmosférica neta de onda larga se 
calculó considerando que la radiación reflejada 
es igual a un 3%, según: 


4 
Hyy = 0.970 (T.,,. +273.16) €, (9) 


donde o = 5.67:10% W m? K*es la constante de 
Stefan-Boltzmamn; T... [*C] es la temperatura 
del aire y Eire eS la emisividad del aire, que se 
estimó según Swinbank (1963), modificado por 
Wunderlich (1968) como: 


Eire = 9 (1+0.17Cb)(T, +273.16) 


alre alre (10) 
donde a, es una constante de proporcionalidad 
igual a 0.937-107. 


Radiación de onda larga emitida por el agua, 
A, 

La radiación de onda larga emitida por el agua 
sigue la ley de Stefan-Boltzmamn: 


H, =t0 (T, +273.16) (11) 
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donde e = 0.97 es la emisividad del agua y T, 
[*C] es la temperatura de la superficie del agua 
y que en el presente estudio se consideró igual 
a la temperatura del agua T, ya que el río se 
encuentra bien mezclado (Link ef al., 2009). 


Tasa de calor por evaporación, H, 


La tasa de calor por evaporación se calcula 
como: 

H, =P, L,E (12) 
donde p, [kg m*] es la densidad del agua, 
L., U kg] es el calor latente del agua igual 
a 1 000 (4 299-2.36T) y E [m s”] es la tasa de 
evaporación según: 

E= Eso, (e, -e,) (13) 
donde F,,., 
viento, 1, [m s”], y toma la forma: 


es una función de la velocidad del 


Ets) =0, +b,w, (14) 
a, [mb ms”] y b, [mb”] son coeficientes sitio- 
específicos de la estación meteorológica. La 
presión de vapor saturado a la temperatura de 
la superficie del agua e, [mb] se calculó con la 
ecuación de Magnus-Tetens: 


IS 0.7858|in(10) 
T,+237.3 
e =e 


Ss 


(15) 


y e, [mb] es la presión de vapor a la temperatura 
del aire que se obtiene a partir de la humedad 
relativa medida de acuerdo con: 


e _HR% 


(16) 
“ * 100 


Tasa de transferencia de calor por conducción, 
H 


C 


La tasa de transferencia de calor por conducción 
se calculó como: 


ad - (TT) 


Ho =P. 1w — Bowen (105) 12) aire (17) 


donde Cu... = 0.61 [mb *C”] es el coeficiente 
de Bowen, P, [mb] es la presión atmosférica, y 
P [mb] es una presión de referencia a nivel del 
mar. La razón entre las presiones se consideró 
igual a la unidad. 

Nótese que los flujos de calor en la 
superficie del agua dependen de cuatro 
variables fundamentales, que son: 1) radiación 
solar incidente, 2) temperatura ambiente, 3) 
velocidad del viento, y 4) humedad relativa del 
aire. 


Formulación discreta del problema 


Las ecuaciones gobernantes se llevaron a la 
forma discreta, aplicando el método de los 
volúmenes finitos. Las derivadas con respecto 
al tiempo se aproximaron utilizando el esquema 
explícito de Euler: 


oT apesar -T! 
dE ON 


(18) 


donde T'es la variable de interés, por ejemplo, 
la temperatura del agua, y t es el tiempo. 
Integrando la ecuación (1) en un volumen de 
control diferencial: 


oT (QT) 0 OT 
ya oa > 09 ED, ON 


(19) 


9T  po(QT T 
tt J As + J D¡AT—nmy9S 


dd 
OX 


(20) 


(29 _ (QT), y, + D¡A 2) 
entra 


Q1) 
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donde n, es un vector normal a la superficie 
diferencial 05 y h es la profundidad del 
escurrimiento en m. 


Transporte advectivo 
Los flujos por advección en un volumen “7” 
se calcularon en las caras “i +?*/,”, siguiendo 
el esquema mostrado en la figura 1. Para su 
cálculo se generalizó la estrategia ultimate 
(Leonard, 1991) a fin de poder utilizar 
parámetros variables en el tiempo y espacio. 
Para el caso en que el intervalo de 
discretización sea constante, y la velocidad 
en la cara del volumen 1,, 1, > 0, se designan 
vecinos de acuerdo con la figura 1 y se verifica 
la monotonicidad, asignando al parámetro 
MON el valor pertinente. A continuación 
se presentan los seis pasos del algoritmo 
propuesto para determinar la temperatura en 
la cara derecha (aguas abajo) del volumen, T,,: 


Paso 1: verificar si la temperatura aumenta 
monotónicamente en el tramo, de modo que se 
cumple: 


Tí <T¿ <Tr <Trg, entonces asigna 


MON =1 (22) 

Si, en cambio, la temperatura disminuye 
monotónicamente en el tramo, de modo que se 
cumple: 


Figura 1. Esquema de discretización espacial en volúmenes 
finitos. 


Tf" > TE > TR > Trg - entonces asigna 


MON=-1 (23) 
Si [MON] +1, se impone T,, = T, y se va al 
paso 6. 


Paso 2: estimar el valor de T,, usando algún 
esquema de alto orden, en este caso quick 
(Leonard, 1979): 


1 6 3 


Toa Ta TE 24 
FR 8 E 8 € 8 R ( ) 


y verificar si se cumple T¿ < T,, < Tr, cuando T 


es monotónicamente creciente o Ti > T,,> Tr, 
cuando T es monotónicamente decreciente. Si 
no se cumple la condición correspondiente, se 


impone T,,= T, y se va al paso 6. 


Paso 3: calcular los valores de referencia: 


G TEST" 
Too = A UG EOL (25) 
REF1 E L a 
C qe Te 
Ther2 = e q > £ (26) 
r r 


donde C, C,y C,, corresponden al número de 
Courant-Friedrichs-Levy (CFL) en las caras 
izquierda (1), derecha (r) y derecha del volumen 
vecino (rr), respectivamente. 


Paso 4: si MON > 0 se asigna: 


Ther = mín (Tagri Trer2) (27) 


y el valor de T,, es limitado por T. abajo, y por 
el menor entre T,,, y T, arriba. 


Paso 5: si MON < 0 se asigna: 


Traer = máx(Tegry , Truro) (28) 
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y el valor de T,, es limitado por T, arriba y por 


el mayor entre T,,, 


y T, abajo. 
Paso 6: actualizar el valor de T de acuerdo con: 


TA TP 4CT, -C Te (29) 


Transporte difusivo 


Los flujos debidos a la difusión que entran y 


MAI 


salen de un volumen “1” se calculan a partir 
del valor del gradiente de la propiedad en sus 
caras “1 + ES Éste se determina a partir de los 
valores en el centro de los volúmenes contiguos 


a dicha cara, ¡e ¡+ 1. Para una discretización 


uniforme: 
9T A 
OX ajo Ax 
El flujo difusivo que entra por 
la cara i - *Y, cuando T,, > T, es 
oT oT 
D, A =- D..,B.,h.,, donde los 
L Ox ÓX ¡1/2 Li-1/21-1/2**i-1/2 


valores del coeficiente de dispersión longi- 
tudinal, ancho de la sección y profundidad en 
la cara se calculan a partir de una interpolación 
lineal de los valores de estas variables en el 
centro de los volúmenes contiguos, ¡e i—1. 


Fuente/sumidero de calor 


El término fuente se considera constante en un 
volumen de control, de manera que: 


e A. (30) 
hC  MCEsbo 


wPw 


con C,, y P, constantes durante el periodo de 
simulación, debido a que los cambios son muy 
pequeños dentro del rango de variación de 
temperatura del río. 


Verificación del modelo 


Advección pura de una propiedad escalar 
conservativa 


Se calculó la advección pura de una campana 
con forma sinusoidal, similar a la que presenta 
el ciclo diario de temperatura en un río. La 


condición inicial es: 


1 A 
TX 
5W, 


donde W, es la longitud de la onda y M 


x>W, 


Tiz,0) = M+ se 


MN , O<x<W, (31) 


la temperatura media. Para efectos de 
comparación, se normalizó con respecto a M. 
En cualquier instante de tiempo, la coordenada 
del pico de concentración, x, e y,, están dadas 
por: 
X¿= Xq +U gd E (32) 
Y.=VY0 


con x, = 0, y, = 1 y t el tiempo. Los intervalos 
de discretización espacial y temporal fueron 
Ax = 2.5 m y At = 0.25 s, respectivamente. La 
longitud de onda, W, = 200 m, y la velocidad 
media de escurrimiento, u,,, = 1 ms”. La figura 
2 muestra la condición inicial y los resultados 
exacto y calculado mediante los esquemas 
upwind, diferencias centradas, quick y ultimate, 
después de 400 s. 

Se observa que los esquemas upwind 
y diferencias centradas presentan una 
difusión numérica significativa. El esquema 
quick conserva de mejor forma el pico de 
concentraciones, pero presenta oscilaciones 
que se propagan en el dominio. La estrategia 
ultimate modificada presenta una solución 
con difusión numérica despreciable y libre de 
oscilaciones. El cuadro 1 muestra la calidad 
de las soluciones numéricas, caracterizada 
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1:5 r r 7 1 1.5 
pe Condición inicial 
..... Solución analítica 
1 a Solución numérica | 1 
2. “Upwind” 


C/Cref 


C/Cref 


Condición inicial 
Solución analítica 
Solución numérica 
“Diferencias centradas] 


1 000 


1,5 , , Ñ , 1,5 , , Ñ , 
0 200 400 600 800 1000 0 200 400 600 800 
Distancia (m) Distancia (m) 
12 1.5 
24385 Condición inicial ===": Condición inicial 
mnmn. Solución analítica =====Solución analítica 
1< Solución numérica lb. se Solución numérica |. 
PA A “Quick” PA $ A 2 “Quick € ultimate” 
O: na 9 
: q a 3 : . $ 3 ] 
si e 4 dE 4 
opos E: O bo? 
9 A b ns Bu $ 4. 
O o E: y ed ys 
pod d Ñ a 9 A 
-0.5+ 4 : -0,5 4 E 
0 $ 3 P 
-1+ -1t : 4 
1,5 , , , , 15 , , , , 
0 200 400 600 800 1 000 0 200 400 600 800 


Distancia (m) 


Distancia (m) 


Figura 2. Advección pura. Condición inicial, y soluciones exacta y numérica para t =400 s. 


1000 


mediante los estimadores: coeficiente de de- 
terminación R?, valor absoluto del promedio 


de las diferencias MAE, error cuadrático 


medio RMS, y eficiencia de Nash y Sutcliffe. 
Las abreviaciones MAE y RMS corresponden 


Cuadro 1. Calidad de las soluciones numéricas para el cálculo de advección pura. 


a las siglas en inglés de los indicadores 
mencionados. 


Se observa que la calidad del esquema 
ultimate es significativamente mejor que la de 
los otros esquemas implementados. 


Esquema de discretización espacial 


Upwind 


Dif. centradas Quick Ultimate 
R? 0.87 0.69 0.99 0.99 
MAE 0.04 0.07 0.01 0.01 
RMS 0.09 0.14 0.03 0.03 
Eficiencia 0.87 0.69 0.99 0.99 
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Difusión pura de un frente de temperatura 


Se calculó la difusión pura de un frente de 
temperatura en un dominio con profundidad 
constante. La condición inicial es: 


0 ] 
E x<X, (83) 
, a , X=X, 


To. es una temperatura de referencia y x, es 
el punto del dominio donde se localiza el frente 
de temperatura. La condición en el borde de 
aguas abajo x = x, es: 


=T (34) 


La solución analítica es (ver e.g., Socolofsky 
y Jirka, 2001): 


(35) 


donde D [m?s] es el coeficiente de difusión. 
Los intervalos de discretización espacial 
y temporal fueron Ax = 2.5 m y At = 0.1 s, 
respectivamente. El tirante fue h = 1 m, el 
coeficiente de difusión D = 10 m's* y x, =1 250 
m. La figura 3 muestra la condición inicial y 
la solución exacta y numérica para t = 500 y 
2500 s. 

Se observa que las soluciones exacta y 
numérica son prácticamente coincidentes. 
El cuadro 2 muestra la diferencia porcentual 
promedio entre ambas soluciones para t = 500 
y 2 500 s. 

La diferencia entre las soluciones exacta y 
numérica para el caso de la difusión pura de 


una frente de temperatura es despreciable. 
Sitio de estudio y antecedentes disponibles 


El modelo desarrollado se aplicó para 


calcular la temperatura en un tramo de 32 km 


Les? Condición inicial] 
[l---- Sol. analítica | 34 
o Sol. numérica 
os <— t=25008 
“+= 0.6 
y 
Ex 
e 
E 04 E 
1 =0s 
0.2 : 
<—1=5008 
+ 


0 500 1.000 1500 2000 2500 
Distancia (m) 


Figura 3. Difusión pura. Condición inicial, y soluciones 
exacta y numérica para £ =500 y 2 500 s. 


Cuadro 2. Diferencia porcentual promedio entre las 
soluciones exacta y numérica para t = 500 y 2 500 s. 


Diferencia porcentual promedio 


100* y abs (Cura 2 a) 


t[s] No 


exacto 


500 0.10 
2 500 0.11 


del río Itata bajo, ubicado en Chile Central. 
La cuenca presenta un clima templado cálido 
con lluvias invernales, caracterizándose 
por la existencia de una estación seca y otra 
lluviosa, con al menos dos meses consecutivos 
de déficit hídrico durante el periodo estival. 
La temperatura ambiente mínima y máxima 
media anual registrada es de 6.8 y 19.9 *C, 
respectivamente, y la precipitación media 
anual en la cuenca es de 1 042.5 mm (Link et 
al., 2009). El tramo de estudio se ubica entre la 
confluencia con el río Ñuble y la localidad de 
Coelemu, entre los paralelos 36* 39 25” y 36" 
28' 10”. La elevación de la cuenca aportante al 
punto final del tramo de estudio varía entre 
7 y 3 815 m, con una pendiente media de 
16.49%. 

En el tramo de estudio existen dos 
estaciones fluviométricas, administradas por 
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la Dirección General de Aguas de Chile, que 
cuentan con registros del caudal medio diario 
en la estación Nueva Aldea, ubicada aguas 
arriba de la confluencia con el río Ñuble y 
Coelemu. En la estación Nueva Aldea, el 
caudal medio diario medido es de 129.8 m's” 
y el caudal instantáneo mínimo y máximo con 
10 años de periodo de retorno, de 1.4 y 2 228.1 
més”, respectivamente. En la estación Coelemu, 
el caudal medio diario medido es de 363.6 m*s* 
y el caudal instantáneo mínimo y máximo con 
10 años de periodo de retorno de 2.1 y 6 885.3 
ns”, respectivamente. El único aporte de agua 
significativo entre las dos estaciones proviene 
del río Ñuble, inmediatamente aguas abajo de 
la estación Nueva Aldea. 

Las variables que controlan el régimen 
térmico en el río Itata bajo se midieron en 
tres estaciones meteorológicas: en la estación 
Nueva Aldea, administrada por la compañía 
Celulosa Arauco y que funciona desde el año 
2006; y en las estaciones Batuco y Trehuaco, que 
se instalaron en el marco del presente estudio. 

La temperatura del agua al inicio y final 
del tramo de estudio se monitoreó mediante 


termistores, con una frecuencia de 15 minutos. 
Las mediciones se utilizaron como condición 
de borde y punto de control para el modelo 
numérico, respectivamente. 

La figura 4 muestra la cuenca del río Itata, 
indicando el tramo de estudio, el río Ñuble; 
la ubicación de las estaciones meteorológicas 
Nueva Aldea, Batuco y Trehuaco; de las 
estaciones pluviométricas Itata en Nueva Aldea 
e Itata en Coelemu; y los puntos de control 
de temperatura del agua en Confluencia y 
Coelemu. 

En febrero de 2009 se realizó una batimetría 
del tramo de estudio mediante un GPS 
diferencial marca Leica y una ecosonda marca 
Garmin, modelo GPSmap 178C. La pendiente 
media del tramo es de 0.3/1 000. El ancho 
superficial y la profundidad con el caudal a 


sección llena, Q varían entre 221-529 m, y 


bankfull 
1.7-2.1 m, respectivamente. 


Resultados 


Se presentan mediciones de la temperatura 
del agua al inicio y al final del tramo de 


América del Sur 
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Figura 4. Cuenca del río Itata (izquierda) y tramo de estudio con puntos de monitoreo (derecha). 
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estudio, así como de las variables que 
controlan el régimen térmico para cuatro 
periodos de doce días de duración cada 
uno, correspondientes a escenarios típicos 
de invierno, primavera, verano y otoño, 
registrados del 2 al 13 de julio de 2009, del 22 
de septiembre al 3 de octubre de 2009, del 2 
al 13 de febrero de 2008, y del 2 al 13 de junio 
de 2009, respectivamente. 

La figura 5 muestra la temperatura medida 
al inicio y final del tramo de estudio, así como 
el caudal correspondiente a los escenarios de 
invierno, primavera, verano y otoño. 

Se observa que la temperatura media 
medida en los escenarios de invierno y 
verano es constante, igual a 9.6 y 24.0 *C, 
respectivamente, y no cambia a lo largo del 
tramo. En los escenarios medidos de otoño 
y primavera se observa un leve aumento 
a lo largo del tramo igual a 0.45 y 0.74 *C, 
respectivamente. La amplitud de la oscilación 
térmica diaria medida en el escenario de 
invierno es despreciable, igual a 0.7 *C, y muy 
significativa en el de verano, igual a 8.4 “C. Se 
observa que sólo en el escenario medido de 
verano dicha amplitud aumenta a lo largo del 
tramo de estudio, a razón de 0.175 “C/km. En 
el escenario de invierno, la temperatura en el 
tramo permanece casi invariable, mientras el 
caudal exhibe una crecida importante del río 
con valores entre 250 y 900 m's*. En el escenario 
de verano, la temperatura exhibe un marcado 
ciclo diario, mientras el caudal permanece muy 
bajo, igual a 8.7 més”. 

La figura 6 muestra las variables medidas 
que controlan el régimen térmico. La radiación 
solar incidente exhibe una marcada fase con 
la temperatura ambiente y con la humedad 
relativa. 

Nótese que el peak de temperaturas del 
agua se alcanza entre dos y cuatro horas 
después que el peak de radiación solar 
incidente. En días nublados, la radiación 
solar máxima en invierno y verano es cercana 
a los 200 (ver días del 1 al 5 en invierno) y 600 
Wm” (ver día 5 en verano); mientras que en 
días descubiertos, la radiación solar máxima 


en invierno y verano es cercana a los 400 (ver 
días 6 y 7, y del 9 al 12 en invierno) y 900 Wm”? 
(ver días del 6 al 12 en verano). Se observa 
una mayor correlación entre la temperatura 
del aire y la radiación en la medida en 
que esta última aumenta, ¡.e. durante los 
escenarios de primavera y verano. En los 
escenarios de primavera y verano se observa 
un aumento del viento durante el día, con 
velocidades máximas entre 3.0-4.5 y 4.0-7.0 
ms”, respectivamente. Durante la crecida del 
escenario de invierno, la velocidad máxima 
del viento alcanzó los 7.0 ms”. 

La figura 7 muestra la temperatura del 
agua medida y calculada al final del tramo 
de estudio para los escenarios de invierno, 
primavera, verano y otoño. 

Se observa que el modelo predice ade- 
cuadamente los valores medidos y reproduce 
su tendencia. El cuadro 3 muestra la calidad de 
la solución numérica, caracterizada mediante 
los estimadores R?, valor absoluto del promedio 
de las diferencias MAE, error cuadrático 
medio RMS y eficiencia para cada uno de los 
escenarios simulados. 

Se nota que en los escenarios de invierno y 
primavera, el R*es relativamente bajo, debido a 
que se producen desfases entre la temperatura 
calculada y la medida. Sin embargo, los valores 
de MAE son bajos, ya que la predicción de los 
valores tiene un error bajo. En los escenarios 
de verano y otoño, los valores de R?son muy 
buenos, ya que el modelo sigue la fase de la 
temperatura medida. El valor de MAE es algo 
mayor que en los otros escenarios analizados, 
pero sigue indicando una calidad buena del 
modelo. Dada la buena concordancia entre los 
resultados de las mediciones y los cálculos, se 
concluye que el modelo propuesto constituye 
una herramienta de simulación con capacidad 
predictiva, que permite la cuantificación de la 
temperatura del agua de un río. 


Conclusión 


Se presentó el desarrollo, verificación y 
aplicación de un modelo numérico para el 
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Figura 5. Temperatura medida al inicio y final del tramo de estudio, y caudal en invierno, primavera, verano y otoño. 
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Figura 7. Temperatura medida y calculada al final del tramo en invierno (arriba, izquierda), primavera 


(abajo, izquierda), verano (arriba, derecha) y otoño (abajo, derecha). 


Cuadro 3. Calidad de la solución numérica respecto de las mediciones al final del tramo en invierno, primavera, verano y otoño. 


Invierno Primavera Verano Otoño 

R? 0.62 0.50 0.81 0.75 
MAE 0.40 0.70 1.09 0.61 
RMS 0.50 0.82 1.33 0.72 
Eficiencia 0.52 0.38 0.78 0.38 


cálculo de la distribución espacio-temporal de 
la temperatura de un río. Se aplicó el método 
de los volúmenes finitos y se generalizó la 
estrategia ultimate para la solución de la ecua- 
ción advección-difusión-reacción, con paráme- 
tros variables que permitieron obtener resul- 
tados libres de difusión numérica. 


El modelo desarrollado se verificó mediante 
el cálculo de dos problemas con solución 
analítica conocida, mostrando su aptitud para 
simular el transporte advectivo y difusivo de 
propiedades escalares en el agua. 

Se midió y calculó la temperatura de un 
tramo de 32 km del río Itata, en cuatro periodos 
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de doce días de duración, correspondientes a 
escenarios de invierno, primavera, verano y 
otoño. 

El modelo desarrollado fue capaz de 
predecir de forma adecuada el complejo com- 
portamiento que exhibió 
medida durante los cuatro escenarios ana- 


la temperatura 


lizados, constituyendo una herramienta de 
simulación con capacidad predictiva. 
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Abstract 


MONSALVE, A., LINK, O. és STEHR, A. The thermal regime of rivers: development, 
verification and application of a numerical model. Water Technology and Sciences (in 
Spanish). Vol. IL, No. 4, October-December, 2012, pp. 41-56. 


The finite volumes method is applied and the ultimate strategy is generalized to solve the 
advection-diffusion-reaction equation with variable parameters, in order to obtain results 
with negligible numeric diffusion. The model is verified through calculations with known 
analytical solutions, showing its capability to simulate diffusive and advective transport of 
scalar properties in water. The model is used to calculate temperature in a 32 km section 
of the lower Itata River, located in Central Chile, measuring the river's temperature and 
variables controlling the thermal regime at several points. Results are presented for four 
12-day periods, corresponding to winter, spring, summer and autumn scenarios. The mean 
temperature measured during the winter and summer scenarios are constant, at 9.6 and 
24.0 *C, respectively, and do not change throughout the section. In the autumn and spring 
scenarios, a slight increase throughout the section is observed, of 0.45 y 0.74 *C, respectively. 
The amplitude of the daily thermal oscillation measured for the winter scenario (0.7 *C) is 
negligible, and very significant (8.4 *C) for the summer scenario. The amplitude is observed 
to increase in the section studied only for the summer scenario, at a rate 0of0.175 “C/km. The 
model developed adequately predicts the complex behaviour of the river's temperature for the 
four scenarios analyzed, constituting a simulation tool with prediction capability. 


Keywords: thermal regime, numerical methods, discretization schemes, Itata River, Chile. 
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